R Notebook

#setwd("C:/Users/erussek/forage_jsp/analysis")
setwd("/Users/evanrussek/forage_jsp/analysis")
rm(list = ls())
library(tidyr)

Attaching package: ‘tidyr’

The following object is masked from ‘package:Matrix’:

    expand

library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(zoo)

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

library(ggplot2)
library(ggpubr)

Loading required package: magrittr

Attaching package: ‘magrittr’

The following object is masked from ‘package:tidyr’:

    extract

library(RcppRoll)
library(knitr)

data <- read.csv('data/run4_data.csv')
head(data)

data <- data %>% ungroup() %>% bind_cols(s_num = group_indices(., subjectID))
data <- data %>% mutate(subj = factor(s_num))

travel_keys = data %>% ungroup() %>% select(travel_key) %>% unique()
clean_subj_data <- function(subj_data, travel_keys){
  travel_key_hard = travel_keys$travel_key[1]
  travel_key_easy = travel_keys$travel_key[2]
  
  subj_data <- data %>%select(round, phase, reward_obs, reward_true, lag, exit, start_reward, n_travel_steps,
                                    travel_key, subjectID, trial_num, s_num, correct_key) %>% group_by(trial_num) %>%
    mutate(press_num = 1:n(), round = as.factor(round), 
           rep_number = case_when(trial_num < 7 ~ 1, TRUE ~ 2)) %>%
    mutate(phase = replace(phase, phase == "Harvest", "HARVEST"))

  subj_data <- subj_data %>% group_by(trial_num) %>% 
    mutate(travel_key = replace(travel_key, travel_key == "", first(travel_key[travel_key != ""]))) %>%
    mutate(travel_key_cond = case_when(travel_key == travel_key_hard   ~ "HARD",
                                       travel_key == travel_key_easy  ~ "EASY")) %>%
    filter(!is.na(travel_key_cond))
}

# clean all the data and bind
datalist <- list()
for (i in 1:10){
  subj_data <- data %>% filter(s_num == i)
  subj_data <- clean_subj_data(subj_data, travel_keys)
  subj_data <- ungroup(subj_data)
  datalist[[i]] <- subj_data
}

Unequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vector

#cdata <- do.call(rbind,datalist)
cdata <- bind_rows(datalist)

# make a function to plot reward on each game... 
plot_subj_reward_v_press <- function(subj_data){
  
  s_num <- subj_data$s_num[1]
  
  # reward plots which show the threshold... 
  rew_plot <- ggplot(subj_data, aes(x = press_num, y = reward_obs, group = round))  +
    #geom_rect(data = t1_data, aes(xmin = press_num - 0.5, xmax = press_num + 0.5, ymin = -Inf, ymax = Inf, fill = round)) +
    geom_point(aes(color = phase)) + facet_grid(n_travel_steps ~ travel_key_cond) + theme(legend.position = "none") + ylim(0,110) + ggtitle(paste("subj: ", s_num)) + ylab('reward collected') + xlab('button press number') + theme_minimal()
  
 # plot_an <- annotate_figure(plot, top = paste('subj: ', s_num))
  plot(rew_plot)
                             
}
for (s in 1:22){
  subj_data <- cdata %>% filter(s_num == s)
  subj_data$round_num <- as.integer(as.character(subj_data$round))

  #subj_data <- subj_data %>% group_by(travel_key_cond, n_travel_steps) %>% mutate(press_num = 1:n()) %>% ungroup()
  plot_subj_reward_v_press(subj_data)  
}

subj_data$round_num 

   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  [41] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  [81] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [121] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [161] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [201] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [241] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [281] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [321] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [361] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [401] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [441] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [481] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [521] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [561] 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [601] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [641] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [681] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [721] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [761] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
 [801] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
 [841] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
 [881] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
 [921] 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
 [961] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 6 6 6
 [ reached getOption("max.print") -- omitted 32240 entries ]

get_trial_exits <- function(thdata){
    # get harvest
  last_phase <- last(thdata$phase)
  
  # go through data.. if last phase was harvest, remove that round... 
  if (last_phase == "HARVEST"){
    last_round <- last(thdata$round)
    # get the last reward ops...
    last_reward_ob <- last(thdata$reward_obs[!is.na(thdata$reward_obs)])
    
    #if (last_reward_ob > 8){
    thdata <- thdata %>% filter(round != last_round)
    #}
  }
  
  # now select harvest data
  thdata <- thdata %>% filter(phase == "HARVEST", !is.na(reward_obs)) # find out why reward_true has so many .na
  
  ## get either last true reward observed
  return_tbl <- thdata %>% 
                  group_by(round) %>% 
                      summarise(s_num = first(s_num), last_reward = last(reward_obs),
                          trial_num = first(trial_num),
                          n_travel_steps = first(n_travel_steps),
                          travel_key_cond = first(travel_key_cond),
                          rep_number = first(rep_number)) %>% ungroup()
  return(return_tbl)
}

#sdata %>% group_by(trial_num) %>%
#  do(get_trial_exits(.))

exit_data <- cdata %>% group_by(s_num, trial_num) %>%
  do(get_trial_exits(.)) %>% ungroup() %>% mutate(subj = as.factor(s_num))

library(ggpubr)
exit_data <- exit_data %>% mutate(round_num = as.integer(as.character(round)))
for (s in 1:22){
  p <- ggplot(exit_data %>% filter(s_num == s), aes(x = round_num, y = last_reward)) + geom_point() + geom_line() +
  facet_grid(n_travel_steps ~ travel_key_cond) + theme_minimal()  + ylim(0,100) + ylab('last reward collected') + xlab('harvest number')
  
  #plot(p)
  #ggarrange(p)
  a <- annotate_figure(p, right = "# travel steps", top = "travel step difficulty", fig.lab = paste('subj: ' , s))
  plot(a)
}

geom_path: Each group consists of only one observation. Do you need to adjust the group aesthetic?

# do some aggregating over exit thresholds... 
# just get the mean for each subject for each timepoint, collapse accross rep number...
mn_exit <- exit_data %>% group_by(s_num, n_travel_steps, travel_key_cond) %>%
  summarise(rep_exit_thresh = mean(last_reward), trial_num = mean(trial_num)) %>% 
  group_by(s_num, n_travel_steps, travel_key_cond) %>%
  summarise(exit_thresh = mean(rep_exit_thresh), trial_num = mean(trial_num)) %>%
  mutate(subj = as.factor(s_num)) %>% ungroup()

Unequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vector

# now just plot this for each subject -- draw lines to show easy - hard effect
ggplot(mn_exit, aes(x = travel_key_cond, y = exit_thresh)) + geom_dotplot(binaxis = 'y', aes(fill = subj), dotsize = 1.2) + geom_line(aes(group = subj, color = subj), size = 1.2) +
  facet_grid(.~n_travel_steps) + theme_minimal() + labs(y = 'Mean Last Reward Collected', x = 'Travel Sequence', subtitle = '# Travel Presses') #theme(legend.position = "none")
ggsave('reward_press2.png')

Saving 7.29 x 4.51 in image

ggplot(mn_exit, aes(x = factor(n_travel_steps), y = exit_thresh)) +
  geom_dotplot(binaxis = 'y', aes(fill = subj), dotsize = 1.2) + 
  geom_line(aes(group = subj, color = subj), size = 1.2) +
  facet_grid(.~travel_key_cond) + 
  theme_gray() + 
  labs(y = 'Mean Last Reward Collected', x = 'Travel Key Condition', subtitle = '# Travel Presses') 
  #theme(legend.position = "none")
ggsave('reward_press.png')

Saving 7.29 x 4.51 in image

# plot the means for each of these...
gmn_exit <- mn_exit %>% group_by(n_travel_steps, travel_key_cond) %>% summarise(gm_thresh = mean(exit_thresh), gsd_thresh = sd(exit_thresh)/sqrt(n()))

ggplot(gmn_exit, aes(x = factor(n_travel_steps), y = gm_thresh, color = travel_key_cond)) +
  geom_line(aes(group = travel_key_cond), size = 2) +
  geom_errorbar(aes(ymin = gm_thresh - gsd_thresh, ymax = gm_thresh+gsd_thresh), width = .1, size = 2) + 
geom_point(size = 4) + ylab('group mean exit threshold') + xlab('# travel steps') + labs(color = 'travel effort') + theme_minimal()
ggsave('group mean exit thresh.png')

Saving 7.29 x 4.51 in image

library(optimx)
library(lmerTest)
exit_model <- lmer(last_reward ~ n_travel_steps + travel_key_cond  + (1  + n_travel_steps + travel_key_cond |subj), data = exit_data,  control = lmerControl(optimizer ="optimx",optCtrl=list(method='nlminb')))
summary(exit_model)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: last_reward ~ n_travel_steps + travel_key_cond + (1 + n_travel_steps +      travel_key_cond | subj)
   Data: exit_data
Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))

REML criterion at convergence: 6956.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.1827 -0.5108  0.0002  0.5046  5.3563 

Random effects:
 Groups   Name                Variance  Std.Dev. Corr       
 subj     (Intercept)         290.72330 17.0506             
          n_travel_steps        0.06055  0.2461  -0.27      
          travel_key_condHARD  68.01256  8.2470  -0.25  0.26
 Residual                      87.82008  9.3712             
Number of obs: 928, groups:  subj, 22

Fixed effects:
                    Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)         54.95291    3.71001 20.77742  14.812 1.63e-12 ***
n_travel_steps      -0.24696    0.06318 19.11207  -3.909 0.000935 ***
travel_key_condHARD -5.06293    1.89212 17.91918  -2.676 0.015465 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) n_trv_
n_trvl_stps -0.309       
trvl_k_HARD -0.256  0.206

library(broom.mixed)
td <- tidy(exit_model, conf.int = TRUE)
tdsome <- td %>% filter(term == 'n_travel_steps' | term == 'travel_key_condHARD')

tdsome <- tdsome %>% mutate(parameter = case_when(term == 'n_travel_steps' ~ 'travel step', term == 'travel_key_condHARD' ~ 'travel effort cond'))

ggplot(tdsome, aes(estimate,parameter,color = parameter)) + 
  geom_point(size = 2)  +
geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), size = 1.2) + theme(legend.position = "none") + theme_minimal() + ylab('fixed effect') + xlab('estimate (# points)') + ggtitle('linear mixed effects model - effect on exit threshold') + theme(axis.text.y = element_text(angle = 90, hjust = .5))
ggsave('exit model.png')

Saving 7.29 x 4.51 in image

# reaction times... 


## a bit more sensible... 
lag_data <- cdata %>% select(s_num, rep_number, travel_key_cond, n_travel_steps, phase, trial_num, lag, correct_key, round, trial_num) %>% mutate(log_lag = log(lag)) %>% group_by(s_num, trial_num, round, phase) %>% slice(-1) %>% ungroup() %>% mutate(subj = as.factor(s_num))

ggplot(lag_data, aes(x=lag)) + geom_histogram(binwidth = 10) + xlim(50,300) + facet_wrap( ~ s_num)

#ggplot(lag_data, aes(x=log_lag))  + geom_histogram(binwidth = .01) + facet_wrap( ~ s_num)

filt_lag <- lag_data %>% group_by(s_num) %>% filter(log_lag < (mean(log_lag) + 2*sd(log_lag)) , log_lag > (mean(log_lag) - 2*sd(log_lag))) %>% ungroup()

ggplot(filt_lag, aes(x=lag)) + geom_histogram(binwidth = 10)  + facet_wrap(~s_num) + xlim(50,350)

ggplot(filt_lag, aes(x=log_lag)) + geom_histogram(binwidth = .02) + facet_wrap(~s_num)

## with a run, plot the mean lag 
round_filt_lag <- filt_lag %>% 
  group_by(s_num, n_travel_steps, travel_key_cond, round, phase) %>%
  summarise(trial_num = first(trial_num), mean_lag = mean(lag), mean_log_lag = mean(log_lag)) %>% mutate(round_num = as.integer(as.character(round)))

for (s in 1:22){
  
  # plot mean lag as a function of round number
  p <- ggplot(round_filt_lag %>% filter(s_num == s), aes(x = round_num, y = mean_log_lag, color = phase)) + geom_point(size = 2) +facet_grid(n_travel_steps ~ travel_key_cond) + geom_line(size = 1.2)  + xlab('harvest/travel number') + ylab('mean log lag') + ggtitle(paste('subj: ', s)) + theme_minimal()
  plot(p)
  
}

NA
NA

library(plotrix)
cond_round_filt_lag <- round_filt_lag %>% 
  group_by(n_travel_steps, travel_key_cond, round_num, phase) %>%
  summarise(mean_log_lag = mean(mean_log_lag), sd_lag = std.error(mean_log_lag))

p <- ggplot(cond_round_filt_lag, aes(x = round_num, y = mean_log_lag, color = phase)) + geom_point(size = 2) +facet_grid(n_travel_steps ~ travel_key_cond) + geom_line(size = 1.2) + geom_linerange(aes(ymin = mean_log_lag - sd_lag, ymax = mean_log_lag + sd_lag)) +
  xlab('harvest/travel number') + ylab('mean log lag') + ggtitle('group: ') +  theme_minimal()

plot(p)

cond_round_filt_lag

trial_filt_lag <- round_filt_lag %>% 
  group_by(s_num, n_travel_steps, travel_key_cond, phase) %>%
  summarise(trial_num = first(trial_num), 
            mean_lag = mean(mean_lag), 
            mean_log_lag = mean(mean_log_lag)) %>% ungroup() %>% mutate(s_num = as.factor(s_num))

ggplot(trial_filt_lag, 
       aes(x = factor(n_travel_steps), y = mean_log_lag, fill = s_num)) + 
  geom_point(aes(color = s_num), size = 2)+ 
  geom_line(aes(group = s_num, color = s_num), size = 1.2) +
  facet_grid(phase ~ travel_key_cond) + theme_minimal() + xlab('# travel steps') + ylab('log lag time')

ggplot(trial_filt_lag, 
       aes(x = travel_key_cond, y = mean_lag, fill = s_num)) + 
  geom_point(aes(color = s_num), size = 2)+ 
  geom_line(aes(group = s_num, color = s_num), size = 1.2) +
  facet_grid(phase ~ n_travel_steps) + theme_minimal()

group_lag <- trial_filt_lag %>% group_by(n_travel_steps, travel_key_cond, phase) %>% summarise(gm_lag = mean(mean_lag), gml_lag = mean(mean_log_lag)) %>% ungroup()

ggplot(group_lag, aes(x = factor(n_travel_steps), y = gml_lag, color = travel_key_cond)) +
  geom_line(aes(group = travel_key_cond), size = 2) + facet_wrap(~phase) + geom_point(size = 2) +xlab('# travel steps') +
  ylab('group mean log lag') + theme_minimal()

NA
NA
NA

harvest_lag <- round_filt_lag %>% ungroup() %>% filter(phase == "HARVEST") %>% mutate(subj = as.factor(s_num))
head(harvest_lag)

hl_model <- lmer(scale(mean_log_lag) ~  n_travel_steps + travel_key_cond + round_num  + trial_num + (1 + n_travel_steps+ travel_key_cond + round_num + trial_num |subj), data = harvest_lag, control = lmerControl(optimizer ="optimx",optCtrl=list(method='nlminb')))

summary(hl_model)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(mean_log_lag) ~ n_travel_steps + travel_key_cond + round_num +  
    trial_num + (1 + n_travel_steps + travel_key_cond + round_num +  
    trial_num | subj)
   Data: harvest_lag
Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))

REML criterion at convergence: 1235.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.5258 -0.4479 -0.0009  0.4449  6.6054 

Random effects:
 Groups   Name                Variance  Std.Dev. Corr                   
 subj     (Intercept)         0.8066735 0.89815                         
          n_travel_steps      0.0003903 0.01976   0.20                  
          travel_key_condHARD 0.1625154 0.40313  -0.02 -0.36            
          round_num           0.0002204 0.01485  -0.19  0.30  0.40      
          trial_num           0.0057660 0.07593  -0.58 -0.60  0.02 -0.14
 Residual                     0.1498717 0.38713                         
Number of obs: 1030, groups:  subj, 22

Fixed effects:
                     Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)         -0.016907   0.197840 21.640619  -0.085 0.932681    
n_travel_steps      -0.002604   0.004533 19.780497  -0.575 0.572061    
travel_key_condHARD  0.083340   0.091057 19.493395   0.915 0.371243    
round_num            0.025677   0.005413 10.465223   4.744 0.000695 ***
trial_num           -0.006761   0.019242 18.808325  -0.351 0.729213    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) n_trv_ t__HAR rnd_nm
n_trvl_stps  0.137                     
trvl_k_HARD -0.040 -0.298              
round_num   -0.215  0.237  0.268       
trial_num   -0.542 -0.533 -0.006 -0.082

library(broom.mixed)
td <- tidy(hl_model, conf.int = TRUE)
tdsome <- td %>% filter(term == 'n_travel_steps' | term == 'travel_key_condHARD' | term == 'round_num' | term == 'trial_num')

tdsome <- tdsome %>% mutate(parameter = case_when(term == 'n_travel_steps' ~ 'travel step', term == 'travel_key_condHARD' ~ 'travel effort cond', term == 'round_num' ~ 'round number', term == 'trial_num' ~ 'trial number'))

ggplot(tdsome, aes(estimate,parameter,color = parameter)) + 
  geom_point(size = 2)  +
geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), size = 1.2) + theme(legend.position = "none") + theme_minimal() + ylab('fixed effect') + xlab('estimate (sd log rt)') + ggtitle('linear mixed effects model - effect on lag - harvest presses') + theme(axis.text.y = element_text(angle = 90, hjust = .5))

travel_lag <- round_filt_lag %>% ungroup() %>% filter(phase == "TRAVEL") %>% mutate(subj = as.factor(s_num))

tl_model <- lmer(scale(mean_log_lag) ~  n_travel_steps*travel_key_cond + round_num   + trial_num + (1 + n_travel_steps*travel_key_cond + round_num + trial_num |subj), data = travel_lag, control = lmerControl(optimizer ="optimx",optCtrl=list(method='nlminb')))

convergence code 1 from optimxboundary (singular) fit: see ?isSingular

summary(tl_model)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(mean_log_lag) ~ n_travel_steps * travel_key_cond + round_num +  
    trial_num + (1 + n_travel_steps * travel_key_cond + round_num +      trial_num | subj)
   Data: travel_lag
Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))

REML criterion at convergence: 919.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.8972 -0.5267 -0.0053  0.5193  4.9564 

Random effects:
 Groups   Name                               Variance  Std.Dev. Corr                         
 subj     (Intercept)                        0.2709010 0.52048                               
          n_travel_steps                     0.0004442 0.02108   0.40                        
          travel_key_condHARD                0.2821369 0.53117   0.28  0.13                  
          round_num                          0.0005459 0.02336  -0.16 -0.17  0.11            
          trial_num                          0.0093260 0.09657  -0.59 -0.37 -0.37  0.15      
          n_travel_steps:travel_key_condHARD 0.0001886 0.01373  -0.14 -0.86 -0.29 -0.28  0.19
 Residual                                    0.1025152 0.32018                               
Number of obs: 1065, groups:  subj, 22

Fixed effects:
                                    Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)                        -0.453891   0.120232 21.484655  -3.775  0.00108 ** 
n_travel_steps                     -0.002340   0.004878 17.764615  -0.480  0.63723    
travel_key_condHARD                 1.248245   0.124035 20.361977  10.064 2.37e-09 ***
round_num                           0.018594   0.006411 11.569233   2.901  0.01376 *  
trial_num                          -0.040759   0.022544 15.132220  -1.808  0.09052 .  
n_travel_steps:travel_key_condHARD  0.010102   0.003814 22.209233   2.649  0.01459 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) n_trv_ t__HAR rnd_nm trl_nm
n_trvl_stps  0.250                            
trvl_k_HARD  0.148  0.205                     
round_num   -0.219 -0.088  0.108              
trial_num   -0.561 -0.359 -0.323  0.086       
n_t_:__HARD  0.023 -0.768 -0.427 -0.198  0.133
convergence code: 1
boundary (singular) fit: see ?isSingular

travel_hard_lag <- round_filt_lag %>% ungroup() %>% filter(phase == "TRAVEL", travel_key_cond == "HARD") %>% mutate(subj = as.factor(s_num))

thl_model <- lmer(scale(mean_log_lag) ~  n_travel_steps + round_num   + trial_num + (1 + n_travel_steps + round_num + trial_num |subj), data = travel_hard_lag, control = lmerControl(optimizer ="optimx",optCtrl=list(method='nlminb')))

summary(thl_model)

Linear mixed model fit by REML. t-tests use Satterthwaite's
  method [lmerModLmerTest]
Formula: 
scale(mean_log_lag) ~ n_travel_steps + round_num + trial_num +  
    (1 + n_travel_steps + round_num + trial_num | subj)
   Data: travel_hard_lag
Control: 
lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))

REML criterion at convergence: 639.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-6.1332 -0.5312 -0.0211  0.4813  5.1650 

Random effects:
 Groups   Name           Variance  Std.Dev. Corr             
 subj     (Intercept)    1.0618889 1.03048                   
          n_travel_steps 0.0001318 0.01148   0.07            
          round_num      0.0024235 0.04923  -0.22  0.39      
          trial_num      0.0126185 0.11233  -0.53 -0.16  0.50
 Residual                0.1582722 0.39783                   
Number of obs: 461, groups:  subj, 22

Fixed effects:
                Estimate Std. Error        df t value Pr(>|t|)
(Intercept)     0.102854   0.233386 20.044025   0.441   0.6641
n_travel_steps  0.012419   0.003756 13.541207   3.307   0.0054
round_num       0.008826   0.013920 13.541148   0.634   0.5366
trial_num      -0.060142   0.030116 11.392098  -1.997   0.0703
                 
(Intercept)      
n_travel_steps **
round_num        
trial_num      . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) n_trv_ rnd_nm
n_trvl_stps -0.055              
round_num   -0.279  0.311       
trial_num   -0.512 -0.271  0.305

library(broom.mixed)
td <- tidy(thl_model, conf.int = TRUE)
tdsome <- td %>% filter(term == 'n_travel_steps' | term == 'travel_key_condHARD' | term == 'round_num' | term == "trial_num")

tdsome <- tdsome %>% mutate(parameter = case_when(term == 'n_travel_steps' ~ 'travel step', term == 'travel_key_condHARD' ~ 'travel effort cond', term == 'round_num' ~ 'round number', term == 'trial_num' ~ 'trial number'))

ggplot(tdsome, aes(estimate,parameter,color = parameter)) + 
  geom_point(size = 2)  +
geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), size = 1.2) + theme(legend.position = "none") + theme_minimal() + ylab('fixed effect') + xlab('sd log lag') + ggtitle('linear mixed effects model - effect on lag - hard travel presses') + theme(axis.text.y = element_text(angle = 90, hjust = .5))

travel_easy_lag <- round_filt_lag %>% ungroup() %>% filter(phase == "TRAVEL", travel_key_cond == "EASY") %>% mutate(subj = as.factor(s_num))

tel_model <- lmer(scale(mean_log_lag) ~  n_travel_steps + round_num   + trial_num + (1 + n_travel_steps + round_num + trial_num |subj), data = travel_easy_lag, control = lmerControl(optimizer ="optimx",optCtrl=list(method='nlminb')))

summary(tel_model)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(mean_log_lag) ~ n_travel_steps + round_num + trial_num +  
    (1 + n_travel_steps + round_num + trial_num | subj)
   Data: travel_easy_lag
Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))

REML criterion at convergence: 889.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.8124 -0.4839 -0.0207  0.4784  3.8018 

Random effects:
 Groups   Name           Variance Std.Dev. Corr             
 subj     (Intercept)    0.912244 0.95511                   
          n_travel_steps 0.001712 0.04138   0.40            
          round_num      0.001649 0.04061  -0.30 -0.36      
          trial_num      0.046682 0.21606  -0.66 -0.65  0.17
 Residual                0.173302 0.41630                   
Number of obs: 604, groups:  subj, 22

Fixed effects:
                Estimate Std. Error        df t value Pr(>|t|)   
(Intercept)     0.103375   0.223186 16.996980   0.463  0.64911   
n_travel_steps -0.010612   0.009613  9.337431  -1.104  0.29727   
round_num       0.035340   0.011024 15.285312   3.206  0.00578 **
trial_num      -0.009579   0.053193  8.845525  -0.180  0.86115   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) n_trv_ rnd_nm
n_trvl_stps  0.291              
round_num   -0.310 -0.230       
trial_num   -0.636 -0.629  0.110

td <- tidy(tel_model, conf.int = TRUE)
tdsome <- td %>% filter(term == 'n_travel_steps' | term == 'travel_key_condHARD' | term == 'round_num' | term == "trial_num")

tdsome <- tdsome %>% mutate(parameter = case_when(term == 'n_travel_steps' ~ 'travel step', term == 'travel_key_condHARD' ~ 'travel effort cond', term == 'round_num' ~ 'round number', term == 'trial_num' ~ 'trial number'))

ggplot(tdsome, aes(estimate,parameter,color = parameter)) + 
  geom_point(size = 2)  +
geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), size = 1.2) + theme(legend.position = "none") + theme_minimal() + ylab('fixed effect') + xlab('sd log lag') + ggtitle('linear mixed effects model - effect on lag - travel easy presses') + theme(axis.text.y = element_text(angle = 90, hjust = .5))

travel_lag <- round_filt_lag %>% ungroup() %>% filter(phase == "TRAVEL") %>% mutate(subj = as.factor(s_num))

tl_model <- lmer(scale(mean_log_lag) ~  n_travel_steps*travel_key_cond + round_num   + trial_num + (1 + n_travel_steps*travel_key_cond + round_num + trial_num |subj), data = travel_lag, control = lmerControl(optimizer ="optimx",optCtrl=list(method='nlminb')))

convergence code 1 from optimxboundary (singular) fit: see ?isSingular

summary(tl_model)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(mean_log_lag) ~ n_travel_steps * travel_key_cond + round_num +  
    trial_num + (1 + n_travel_steps * travel_key_cond + round_num +  
    trial_num | subj)
   Data: travel_lag
Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))

REML criterion at convergence: 919.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.8972 -0.5267 -0.0053  0.5193  4.9564 

Random effects:
 Groups   Name                               Variance  Std.Dev. Corr                   
 subj     (Intercept)                        0.2709010 0.52048                         
          n_travel_steps                     0.0004442 0.02108   0.40                  
          travel_key_condHARD                0.2821369 0.53117   0.28  0.13            
          round_num                          0.0005459 0.02336  -0.16 -0.17  0.11      
          trial_num                          0.0093260 0.09657  -0.59 -0.37 -0.37  0.15
          n_travel_steps:travel_key_condHARD 0.0001886 0.01373  -0.14 -0.86 -0.29 -0.28
 Residual                                    0.1025152 0.32018                         
      
      
      
      
      
      
  0.19
      
Number of obs: 1065, groups:  subj, 22

Fixed effects:
                                    Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)                        -0.453891   0.120232 21.484655  -3.775  0.00108 ** 
n_travel_steps                     -0.002340   0.004878 17.764615  -0.480  0.63723    
travel_key_condHARD                 1.248245   0.124035 20.361977  10.064 2.37e-09 ***
round_num                           0.018594   0.006411 11.569233   2.901  0.01376 *  
trial_num                          -0.040759   0.022544 15.132220  -1.808  0.09052 .  
n_travel_steps:travel_key_condHARD  0.010102   0.003814 22.209233   2.649  0.01459 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) n_trv_ t__HAR rnd_nm trl_nm
n_trvl_stps  0.250                            
trvl_k_HARD  0.148  0.205                     
round_num   -0.219 -0.088  0.108              
trial_num   -0.561 -0.359 -0.323  0.086       
n_t_:__HARD  0.023 -0.768 -0.427 -0.198  0.133
convergence code: 1
boundary (singular) fit: see ?isSingular

#term == 'travel_key_condHARD'

td <- tidy(tl_model, conf.int = TRUE)
tdsome <- td %>% filter(term == 'n_travel_steps'  | term == 'round_num' | term == "trial_num" | term == "n_travel_steps:travel_key_condHARD")

tdsome <- tdsome %>% mutate(parameter = case_when(term == 'n_travel_steps' ~ '# step', term == 'travel_key_condHARD' ~ 'effort', term == 'round_num' ~ 'round', term == 'trial_num' ~ 'trial', term == 'n_travel_steps:travel_key_condHARD' ~ '#step x effort'))

ggplot(tdsome, aes(estimate,parameter,color = parameter)) + 
  geom_point(size = 2)  +
geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), size = 1.2) + theme(legend.position = "none") + theme_minimal() + ylab('fixed effect') + xlab('sd log lag') + ggtitle('linear mixed effects model - effect on lag - travel') + theme(axis.text.y = element_text(angle = 90, hjust = .5))

td